{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "38c44e25-012b-4631-9c99-48149d5dae46",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a57ae696-0c96-4129-880b-ae02012b3c21",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 光变曲线模拟主函数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "265cca4d-b7c8-4b5c-b4d6-0b2db0047914",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "aea0335a-87e1-4527-8714-eabfb87e54ea",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: AstropyDeprecationWarning: The update_default_config function is deprecated and may be removed in a future version. [sncosmo]\n"
     ]
    }
   ],
   "source": [
    "from my_simu_pkg.lc_simu import lc_simu\n",
    "from my_simu_pkg.lc_simu import load_csst_filter, reset_csst_filter\n",
    "import sncosmo\n",
    "import datetime as dt_\n",
    "import random as rnd\n",
    "source = \"salt2-extended\"\n",
    "model = sncosmo.Model(source=source)\n",
    "csst_bands = load_csst_filter()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "773d428a-c061-4dc1-ad8a-1e2128d301e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "dust = sncosmo.CCM89Dust()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "af7d320b-ba0b-4469-826f-04bc20b5518c",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = sncosmo.Model(source = source, effects=[dust, dust], effect_names=['host', 'mw'], effect_frames=['rest', 'obs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d5f7ee08-c097-418c-b1b9-1f4104e23f97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# class DataQualityError(Exception):\n",
    "#     pass\n",
    "from sncosmo.fitting import DataQualityError\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "56c4a6b5-f951-41ce-ae8e-893dd0c1c0aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import my_simu_pkg as simu\n",
    "import numpy as np\n",
    "import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fc284d0b-7dd2-4ef6-a707-153b72baf5f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = simu.generate_params.get_random_params(t0_min=56150, t0_max=56150+365.25*2, csst_duration = 365.25*2) # 配置CSST巡天起止时间\n",
    "params = np.asarray(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "12f375d5-83ec-4ccc-849d-9e841f7169b6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6665, 7)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8268ae83-86b0-42d4-b63e-931b905352e5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b6cdd8ff-2a0c-45b4-93b2-4bd1e8a75671",
   "metadata": {},
   "outputs": [],
   "source": [
    "names = 'zdist, t0dist, x0dist, x1dist, cdist, Mdist, mbdist'.split(', ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "719f887c-1a19-4cc8-acf7-7b7c9ac74ea7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAHSCAYAAAAT0iZvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABZgElEQVR4nO3de5xVZd3//9dHBEaD5Ji3MtLgHaaIOOYongPNExbYnSnmATzE7SntLktMUzMtLL/5s1vviDwAlYKa5dxpHpGwEnUwRBRNRLwZQsHhoKYo4Of3x7oGN8OemTV79tprH97Px2M/Zu112p917TXX/uxrX+ta5u6IiIiIiEj7tkk7ABERERGRUqHkWUREREQkJiXPIiIiIiIxKXkWEREREYlJybOIiIiISExKnkVEREREYto27QAA+vXr5zU1NWmHISKSk3nz5r3l7v3TjqNQVGeLSCnrbJ1dFMlzTU0NDQ0NaYchIpITM3s97RgKSXW2iJSyztbZ6rYhIiIiIhKTkmcRERERkZiKotuGiORXzcT7t3i+dNJxKUUiIiJJUD2fnqJNnjds2EBjYyPr169PO5TEVVVVUV1dTdeuXdMORURERMpEJeVS2SSVXxVt8tzY2EjPnj2pqanBzNIOJzHuTlNTE42NjQwaNCjtcERERKRMVEoulU2S+VXRJs/r16+viDfbzOjbty+rVq1KOxRJQcuf3aD9n95y2SZf9DOhiEjpqJRcKpsk86uivmCwUt7sSjlOERERKaxKzjGSOvaibXkuBVOnTqWhoYGbbrqJyZMns/3223P66adnXXf27Nl069aNgw46qMBRihQXtV5LuSjmc7mYY5PKsnbtWu644w7OO+88AKZNm8Y111wDwOWXX864ceO22qbY86uSSZ6z/VTdGfmuSM4555w2l8+ePZsePXooeRYREekgfRnIjzRyqbVr1/I///M/nHfeeaxevZof/OAHNDQ0YGbsu+++jB49mt69e7e6fTHmVyWTPKdh8uTJTJ48GYB169ZRU1PD6aefzo9//GN69erF3nvvTffu3QG46qqr6NGjBxdffDE///nPmTx5Mttuuy1Dhgxh0qRJTJ48mS5duvCb3/yG//7v/+bQQw9N89BEtpLLh1NaH2j6IBURKQ0TJ07k1Vdfpba2ln/9618ceeSR9OnTB4AjjzySBx98kJNPPpnbb7+9ZPKrdpNnM7sN+CKw0t2HhnlXAV8Hmnthf8/dHwjLLgXOAjYBF7r7QwnEXRDnnHMO55xzDhs2bODwww/nzDPP5LLLLmPevHnssMMOjBw5kn322Wer7SZNmsRrr71G9+7dWbt2Lb169eKcc87Z/OaLiEikvS9CpXZRbSGU+/FJeZk0aRILFy5k/vz5XH/99VsMm1ddXc3y5ctZsWIFV155ZcnkV3FanqcCNwHTW8y/wd2vz5xhZkOAscCewM7Ao2a2m7tvykOsqbnooos4/PDD6dmzJyNGjKB///4AnHTSSfzjH//Yav1hw4ZxyimncPzxx3P88ccXOFqpRPn+KS6pfYoUC/16IeUmW539p3G7FuS1FzSu3eL5sOpeHdr+qaeeKqn8qt3k2d3nmFlNzP2NAWa4+wfAa2a2GNgfeDL3ENM1depUXn/9dW666Sbq6+tjbXP//fczZ84c/vd//5drr72W559/PuEoRdpW6i1VSuSlLTo/yk+lv6elXmdnJtPLV7y9eXrAgAHMnj178/PGxkZGjBgRe7/Fkl91ps/zBWZ2OtAAfNvd1wADgLkZ6zSGeVsxswnABICBAwd2IozkzJs3j+uvv54nnniCbbbZhuHDh3PRRRfR1NTEJz/5Se6++2723nvvLbb56KOPWLZsGSNHjuSQQw5hxowZvPvuu/Ts2ZO33367lVeSSpHLB0K5f4ioBVBEpH2l+lnwiR49eOeddwA4+uij+d73vseaNWsAePjhh/nxj3/MBx98UFL5Va7J8y+AHwIe/v4/4MyO7MDdpwBTAOrq6jzHOBJ10003sXr1akaOHAlAXV0dV111FQceeCC9evWitrZ2q202bdrEqaeeyrp163B3LrzwQnr16sWXvvQlTjjhBO677z5dMChtKtUKUiQfdP4Xj1Jv/ZTkteyukU2v3n04+OCDGTp0KMceeyzf//732W+//QC44oorNl88GCe/+vKJY3n37bdxd7467uup5Vc5Jc/u/mbztJn9CvhjeLoc2CVj1eowr9PS+Ie9/fbbs84/44wztpp31VVXbZ7+y1/+stXy3XbbjQULFuQtNpFSoERISkUxJ4r6P5J8qb/g4C2ex+mbHCdBbs8dd9yxxfMzzzxz836b/55xxhnt5lfT7n1wq+Vp5Fc5Jc9mtpO7rwhPvwwsDNP1wB1m9jOiCwYHA093OkqRTijmD8Vyk48PeSUKIpWlUP/z+iyQfIkzVN2dwAign5k1AlcCI8yslqjbxlLgPwHc/QUzuwt4EdgInF/qI21IZchHv1slfaVFH6SSb3HqgFI/x0rpGoW06uRSKiPJTZzRNk7OMvvWNta/Fri2M0GJiEjuWhmfvw8wE6ghavQ40d3XmJkBNwKjgPeA8e7+bBpxixQbfcmWbIr6DoPuTlSvlzf3orxeUtqhlubioQ+4rUxl6/H5JwKPufskM5sYnl8CHEvUxW4wMJzogvDhBY1WJCWV8KtjMeRS+eg3nYuk8quiTZ6rqqpoamqib9++qb/pSXJ3mpqaqKqqSjuUslHsFZlI0loZn38MURc8gGnAbKLkeQww3aNPmblm1qvFdS1SBiqhS4lsrVJyqWySzK+KNnmurq6msbGRVatWtb9yiauqqqK6ujrtMKTC6EtGxdkxIyF+A9gxTA8AlmWs1zw+f8Unz+q72j6VUXEbM/UlvjG8N5/u9TrG1snzone2a3cfb655P4nQtpJLLO1tk1R+VbTJc9euXRk0aFDaYYiIlB13dzPr0O+ZpXBjKykNSX1xL+UGgaRif/uDj7h2TlOry+N82Tm2QOWaSyxpfVkr2uRZJBv1bRXJ2ZvN3THMbCdgZZgfa3z+UrixVSUr5cRRpNQoeRYRqQz1wDhgUvh7X8b8C8xsBtGFguvKsb9zsSaX5dYgUG7HU+70fuVGybMUjPrGSZqKNXlKQivj808C7jKzs4DXgRPD6g8QDVO3mGiouq1v8SV5U8znYTHHJoWj86B9Sp5FpCxV8gdAK+PzAxyRZV0Hzk82ovwq9/e23I+v1On9ESXPIjGpwhQR6bhKrzvVNaL8KHkWERGRklLpCbmkS8mzVKT2Kl5VzCLSkuoFEQHYJu0ARERERERKhVqeRURERApIv2JsrZRG5FLLs4iIiIhITGp5lpKnb/AiIiLlpZg/25U8i4iIiFDcCZsUDyXPEksSfZHiVFLF3OdJRDpGiYmIlIN2+zyb2W1mttLMFmbM62Nmj5jZK+Fv7zDfzOznZrbYzBaY2eeSDF5EREREpJDiXDA4FTimxbyJwGPuPhh4LDwHOBYYHB4TgF/kJ0wRERERkfS1mzy7+xxgdYvZY4BpYXoacHzG/OkemQv0MrOd8hSriIiIiEiqcu3zvKO7rwjTbwA7hukBwLKM9RrDvBWI5EB9JEVERKSYdHqcZ3d3wDu6nZlNMLMGM2tYtWpVZ8MQEREREUlcrsnzm83dMcLflWH+cmCXjPWqw7ytuPsUd69z97r+/fvnGIaIiIiISOHk2m2jHhgHTAp/78uYf4GZzQCGA+syundIGcvWvULDzImIiEi5aTd5NrM7gRFAPzNrBK4kSprvMrOzgNeBE8PqDwCjgMXAe8AZCcQsIiIiIpKKdpNndz+5lUVHZFnXgfM7G5SIiIiISDHq9AWDIiIiIiKVQrfnlsRomDkREREpN2p5FhERERGJSS3PklMLsVqVRUREpBKp5VlEREREJCYlzyIiIiIiMSl5FhERERGJScmziIiIiEhMumCwAuliP5HKZWZLgXeATcBGd68zsz7ATKAGWAqc6O5r0opRRKSYqeVZRKTyjHT3WnevC88nAo+5+2DgsfBcRESyUPIsIiJjgGlhehpwfHqhiIgUN3XbqADqpiEiGRx42Mwc+KW7TwF2dPcVYfkbwI6pRSciUuSUPIuIVJZD3H25mX0KeMTMXspc6O4eEustmNkEYALAwIEDCxOpiEgRUrcNEZEK4u7Lw9+VwO+B/YE3zWwngPB3ZZbtprh7nbvX9e/fv5Ahi4gUFSXPIiIVwsw+YWY9m6eBo4CFQD0wLqw2DrgvnQhFRIqfum2IiFSOHYHfmxlE9f8d7v6gmT0D3GVmZwGvAyemGKOISFFT8iwiUiHcfQmwd5b5TcARhY9IRKT0dCp51mD7IiIiIlJJ8tHyPNLd38p43jzY/iQzmxieX5KH15EsWg5Dt3TScSlFIiIiIlL+krhgUIPti4iIiEhZ6mzy3DzY/rwwBihosH0RERERKVOd7baR02D7oAH3RURERKT0dKrlOdfB9sM2GnBfREREREpKzi3PYYD9bdz9nYzB9q/m48H2J6HB9jul5cWA0P4Fgdm2EREREZH86Ey3DQ22LyIiIiIVJefkWYPtp0MtyyIiIiLpSWKoOhERERGRsqTkWUREREQkJiXPIiIiIiIx5eP23BJDnJEz1J9ZREREpLgpeU6RkmURERGR0qJuGyIiIiIiManlOSFqVRYREREpP2p5FhERERGJSS3POVCrsoiIiEhlUsuziIiIiEhMSp5FRERERGJSt40Y1E1DREREREAtzyIiIiIisSl5FhERERGJScmziIiIiEhM6vOchfo4i4iIiEg2ankWEREREYlJybOIiIiISEyJddsws2OAG4EuwC3uPimp1+oMddEQESmdOltEJG2JJM9m1gW4GTgSaASeMbN6d38xidfrCCXLIiJbKuY6W0Sk2CTV8rw/sNjdlwCY2QxgDJDXiliJsIhIXhSkzhYRKQdJJc8DgGUZzxuB4ZkrmNkEYEJ4+q6ZvdzKvvoBb+U9wvwp5vgUW+6KOT7Flrs247Prct7vp3Pesjjks84ulGI91xRXxyiujlFcGdqps9uKqVN1dmpD1bn7FGBKe+uZWYO71xUgpJwUc3yKLXfFHJ9iy12xx1fM4tbZhVKs76Xi6hjF1TGKK74kY0pqtI3lwC4Zz6vDPBERKT6qs0VEYkoqeX4GGGxmg8ysGzAWqE/otUREpHNUZ4uIxJRItw1332hmFwAPEQ17dJu7v5Dj7ormZ8JWFHN8ii13xRyfYstdsceXijzX2YVSrO+l4uoYxdUxiiu+xGIyd09q3yIiIiIiZUV3GBQRERERiUnJs4iIiIhITKkmz2Z2jJm9bGaLzWxiluXdzWxmWP6UmdVkLLs0zH/ZzI5OIbZvmdmLZrbAzB4zs09nLNtkZvPDI+8X3cSIbbyZrcqI4eyMZePM7JXwGJfv2GLGd0NGbP8ws7UZy5Iuu9vMbKWZLWxluZnZz0PsC8zscxnLEi27GLGdEmJ63sz+ZmZ7ZyxbGubPN7OGFGIbYWbrMt67KzKWtXk+FCi+72TEtjCcZ33CskTLTtrWWvmb2TfM7CUze8HMfhLmHWlm88L688zs8Iz19w3zF4f/YQvz+5jZI+H/9hEz651AXPtnnF/PmdmXM9bPev5bdHHmU2H+TIsu1MxrXBnLBprZu2Z2cTHEZWY1ZvZ+RplNzlg/tfcxzB9mZk+G+c+bWVXacVlU98/PeHxkZrVFEFdXM5sW1l9kZpdmrJ+386uDMXUzs9vD+s+Z2YiM9fNXVu6eyoPoopRXgV2BbsBzwJAW65wHTA7TY4GZYXpIWL87MCjsp0uBYxsJbB+mz22OLTx/N+VyGw/clGXbPsCS8Ld3mO5d6PharP8NoouTEi+7sP/DgM8BC1tZPgr4E2DAAcBTBSy79mI7qPk1gWObYwvPlwL9Uiy3EcAfO3s+JBVfi3W/BMwqVNnp0e77sVX5h/r1UaB7eP6p8HcfYOcwPRRYnrHN0+F/1sL/8LFh/k+AiWF6InBdAnFtD2wbpncCVhJdkN/q+Q/cBYwN05OBc/MdV8bye4C7gYvD81TjAmraqEvSfB+3BRYAe4fnfQm5RZpxtVhnL+DVIimvrwEzMv4Hlob3Nq/nVwdjOh+4vXkeMA/YJt9llWbL8+bbwbr7h0Dz7WAzjQGmhel7gCPCN4UxRG/YB+7+GrA47K9gsbn74+7+Xng6l2hc1EKIU26tORp4xN1Xu/sa4BHgmJTjOxm4M88xtMrd5wCr21hlDDDdI3OBXma2EwUou/Zic/e/hdeGwp5zccqtNZ05X2PrYHwFPeckJ+cCk9z9AwB3Xxn+/t3d/xnWeQHYzqJfKHcCPunucz36BJwOHB/Wy/wcmZYxP59xvefuG8M6VUDzlfhZz//wOXY40edaYnEBmNnxwGtE5dUs9biySft9BI4CFrj7c2F+k7tvKoK4Mp1M9H4VQ3k58Akz2xbYDvgQeJvCnF+txTQEmJUxby1Ql++ySjN5znY72AGtrRMqpnVE3wTjbJt0bJnOIvoW06zKzBrMbG6ouPIpbmxfsegn/nvMrPnmB0mXW4dew6KuLoMIJ3qQZNnF0Vr8hSi7jmh5zjnwsEU/ZU9oZZukHRh+JvuTme0Z5hVVuZnZ9kRfen6XMbsYyq6SZSv/3YBDw8+7fzaz/bJs9xXg2fDhOYDo3GqWeZ7t6O4rwvQbwI5JxGVmw83sBeB54JzwmdXa+d8XWJuRcHfk/yJ2XGbWA7gE+EGLfaQaVzDIzP4e5h+aEVea7+NugJvZQ2b2rJl9t0jiynQSH3/5Tzuue4B/ASuA/wOud/fV5P/86khMzwGjzWxbMxsE7Et0A6i8llVqt+cuF2Z2KlAHfD5j9qfdfbmZ7QrMMrPn3f3VAob1v8Cd7v6Bmf0n0Tepw9vZJg1jgXvcfVPGvLTLruiZ2Uii5PmQjNmHhHL7FPCImb0UWmML5Vmi9+5dMxsF/AEYXMDXj+tLwF9DBd8s7bKrdFuVP9FnUx+in1j3A+4ys11DixHhy9l1RC2Fsbm7m5m3v2bH43L3p4A9zWwPYJqZ/an1XXdK7LiAq4Abwv9lQuHkFNcKYKC7N5nZvsAfMr5wtyup9zHMPyTMew94zMzmETXcpRZXxnk/HHjP3bNe21HouIhamDcBOxN1Z3zCzB7tSGwJxHQbsAfQALwO/C3EGEvcskqz5TnO7WA3r2PRzwI7AE0xt006NszsC8BlwOjmnw4A3H15+LsEmE3UT69gsYWfmprjuYXom1esbQsRX4axtPj5POGyi6O1+Ivi9sVmNozoPR3j7k3N8zPKbSXwe/Lbjald7v62u78bph8AuppZP4qk3DK0dc6lUnaVrpXybwTuDUnp08BHQD8AM6sO652e8cV6OVt2Y8o8z9606Cfb5p+52+w6kGtcGdstAt4l9Mkm+/nfRNQlbNsW8/Md13DgJ2a2FPgm8D2LboaTalwedblsCuvPI+ofuxvpv4+NwBx3f8ujbpkPEF1LkXZczVrWX2nH9TXgQXffENb/K1FjYl7Prw6eWxvd/b/cvdbdxwC9gH/ku6zSTJ7j3A62HhgXpk8gusjHw/yxFvV1G0TUwvV0IWMzs32AXxIlzpn9y3qbWfcw3Q84GHixwLHtlPF0NLAoTD8EHBVi7E3UavNQHmOLFV+IcXeib6pPZsxLuuziqAdOt8gBwDqPfs4pRNm1ycwGAvcCp7n7PzLmf8LMejZPh9g61DKRh9j+zWzzlcv7E9UtTRTRbZ/NbAeiX4juy5iXetlVsjbK/w9EFwRhZrsRXXT0lpn1Au4nurjnr837Cf+jb5vZAeE8PJ2P3+fMz5FxGfPzGdeg5oTAou5ouxNd5JT1/A+fY48Tfa4lFpe7H+ruNe5eA/x/wI/c/aa04zKz/mbWJczflegzfEna7yNRnb6XmW0f3s/PAy8WQVyY2TbAiYT+zpD+eU/UVePwjPUPAF4ij+dXDufW9mE9zOxIYKO75/09bPfKyyQfRCMb/IPoW+dlYd7VRAkpRBde3E10QeDTwK4Z214WtnuZcMVkgWN7FHgTmB8e9WH+QUR93p4Lf89KIbYfE10c8hzRibp7xrZnhvJcDJyRxvsanl9F1Nk/c7tClN2dRD8ZbiD65noWcA5RH0WIrsK9OcT+PFBXqLKLEdstwJqMc64hzN81lNlz4X2/LIXYLsg45+YCB7V1PhQ6vrDOeMKV4RnzEi87Pdp837KWP9EH4W+IPiSfBQ4P8y8n6mM5P+PRfJV9XVj/VeAm2HwH3b7AY8ArRPV2nwTiOi2sNz/MPz5jX1nP//AaT4f65G7CqAH5jKvFtlcRRttIOy6i/uqZ5fWljH2l9j6GZaeGdRcCPymiuEYAc7PsK83zvkc4R14gauj6Tr7PrxxiqiHKCxeF4/50EmWl23OLiIiIiMSkOwyKiIiIiMSk5FlEREREJKaiGKquX79+XlNTk3YYIiI5mTdv3lvu3j/tOApFdbaIlLLO1tlFkTzX1NTQ0NDQ/ooiIkXIzF5PO4ZCUp0tIqWss3W2um2IiIiIiMSk5FlEREREJCYlzyIiIiIiMRVFn+dc1Uy8f4vnSycdl1IkIuVtw4YNNDY2sn79+rRDSVVVVRXV1dV07do17VBEEtfyMxa2/pzV53BxUp0dSarOLunkWUQKo7GxkZ49e1JTU0O4C3fFcXeamppobGxk0KBBaYcjUjbiJOnSMaqzk62z1W1DRNq1fv16+vbtW7GVMICZ0bdv34pvyRGR4qc6O9k6W8mziMRSyZVwM5WBiJQK1VfJlYG6bYiIiEhe5NJPWqTUKHkWkQ7L94dfrv0bp02bxjXXXAPA5Zdfzrhx47ZaZ+rUqTQ0NHDTTTcxefJktt9+e04//fSs+5s9ezbdunXjoIMOyikeEZFiVCx19jHHHMPcuXM55JBD+OMf/5h1ndmzZ3P99dfzxz/+kfr6el588UUmTpyYdd358+fzz3/+k1GjRuUUT66UPItISVq9ejU/+MEPaGhowMzYd999GT16NL179251m3POOafNfc6ePZsePXooeRYRScB3vvMd3nvvPX75y1/GWn/06NGMHj261eXz58+noaGh4Mmz+jyLSNF75plnGDZsGOvXr+df//oXe+65JzfffDNHHnkkffr0oXfv3hx55JE8+OCDANx+++3stttu7L///vz1r3/dvJ+rrrqK66+/HoCf//znDBkyhGHDhjF27FiWLl3K5MmTueGGG6itreWJJ55I5VhFilnNxPu3eJT660gystXZCxcu5IgjjqBnz55brf/ggw+y++6787nPfY5777138/ypU6dywQUXAHD33XczdOhQ9t57bw477DA+/PBDrrjiCmbOnEltbS0zZ84s2PGp5VlEit5+++3H6NGjufzyy3n//fc59dRT6dq1K7vsssvmdaqrq1m+fDkrVqzgyiuvZN68eeywww6MHDmSffbZZ6t9Tpo0iddee43u3buzdu1aevXqxTnnnEOPHj24+OKLC3l4eWVmVcAcoDtRHX+Pu19pZoOAGUBfYB5wmrt/aGbdgenAvkATcJK7L00leBEpC9nq7KFDh2Zdd/369Xz9619n1qxZfOYzn+Gkk07Kut7VV1/NQw89xIABA1i7di3dunXj6quv3twtr5DU8iwiJeGKK67gkUceoaGhge9+97utrvfUU08xYsQI+vfvT7du3VqtiIcNG8Ypp5zCb37zG7bdtqzaET4ADnf3vYFa4BgzOwC4DrjB3T8DrAHOCuufBawJ828I64mIdErcOvull15i0KBBDB48GDPj1FNPzbrewQcfzPjx4/nVr37Fpk2bkgo7FiXPIlISmpqaePfdd3nnnXdYv349AwYMYNmyZZuXNzY2MmDAgNj7u//++zn//PN59tln2W+//di4cWMSYRecR94NT7uGhwOHA/eE+dOA48P0mPCcsPwI0xhXkkfqglGZWtbZnTV58mSuueYali1bxr777ktTU1MeosyNkmcRKQn/+Z//yQ9/+ENOOeUULrnkEo4++mgefvhh1qxZw5o1a3j44Yc5+uijGT58OH/+859pampiw4YN3H333Vvt66OPPmLZsmWMHDmS6667jnXr1vHuu+/Ss2dP3nnnnRSOLr/MrIuZzQdWAo8ArwJr3b35G0Ij0PxNYwCwDCAsX0fUtaPlPieYWYOZNaxatSrhIxCRUteyzm7N7rvvztKlS3n11VcBuPPOO7Ou9+qrrzJ8+HCuvvpq+vfvz7Jly1Krs8vqt0oRKYxC3zp3+vTpdO3ala997Wts2rSJgw46iPnz5/P973+f/fbbD4h+IuzTpw8QXRh44IEH0qtXL2pra7fa36ZNmzj11FNZt24d7s6FF15Ir169+NKXvsQJJ5zAfffdx3//939z6KGHFvIw88bdNwG1ZtYL+D2wex72OQWYAlBXV+ed3Z+kr2UrcDHfEruUYi1GxVBnz5o1iyuvvJKXXnqJd999l+rqam699VaOPvpopkyZwnHHHcf222/PoYcemjUh/s53vsMrr7yCu3PEEUew9957M3DgQCZNmkRtbS2XXnppq9308s3c068D6+rqvKGhocPb6Z9JpDAWLVrEHnvskXYYRSFbWZjZPHevSymkNpnZFcD7wCXAv7n7RjM7ELjK3Y82s4fC9JNmti3wBtDf2/hwyLXOluLS3mdoMXex0Od921RnfyyJOjvv3TbMbA8zm2xm95jZufnev4iItM7M+ocWZ8xsO+BIYBHwOHBCWG0ccF+Yrg/PCctntZU4i4hUuljJs5ndZmYrzWxhi/nHmNnLZrbYzCYCuPsidz8HOBE4OP8hi4hIG3YCHjezBcAzwCPu/keiludvmdlioj7Nt4b1bwX6hvnfArLfyktERID4fZ6nAjcRjQUKRBekADcTtWo0As+YWb27v2hmo4FzgV/nN1wRSYu7U+mDMJRCg6y7LwC2Gtja3ZcA+2eZvx74agFCE5ECUp2dXJ0dq+XZ3ecAq1vM3h9Y7O5L3P1DosH3x4T16939WOCUfAYrIumoqqqiqampJJLHpLg7TU1NVFVVpR2KiEibVGcnW2d3ZrSNzcMbBY3AcDMbAfwH0d2tHmhtYzObAEwAGDhwYCfCEJGkVVdX09jYSKUPUVZVVUV1dXXaYYiItEl1diSpOjvvQ9W5+2xgdoz1NOyRSIno2rUrgwYNSjsMERGJQXV2sjoz2sZyYJeM59VhnoiIiIhIWepM8vwMMNjMBplZN2As0ZBHIiIiIiJlKe5QdXcCTwKfNbNGMzsr3Mb1AuAhojFE73L3F5ILVUREREQkXbH6PLv7ya3Mf4A2LgoUERGR0lDMdxQUKSZ5v2BQREREiku5JcbZjke37JZCyfvtuUVEREREypWSZxERERGRmJQ8i4iIiIjEpORZRERERCQmJc8iIiIiIjFptA0REREpOy1H5NBoHJIvSp5FRESk5JXbcHxSvNRtQ0REREQkJiXPIiIiIiIxqduGiEgZMbNdgOnAjoADU9z9RjPrA8wEaoClwInuvsbMDLgRGAW8B4x392fTiF1yp/69IoWj5FlEpLxsBL7t7s+aWU9gnpk9AowHHnP3SWY2EZgIXAIcCwwOj+HAL8JfKWHq/yuSHHXbEBEpI+6+ornl2N3fARYBA4AxwLSw2jTg+DA9BpjukblALzPbqbBRi4iUDiXPIiJlysxqgH2Ap4Ad3X1FWPQGUbcOiBLrZRmbNYZ5IiKShZJnEZEyZGY9gN8B33T3tzOXubsT9YfuyP4mmFmDmTWsWrUqj5GKiJQW9XkWESkzZtaVKHH+rbvfG2a/aWY7ufuK0C1jZZi/HNglY/PqMG8L7j4FmAJQV1fXocRbpBhk6weuCyslF2p5FhEpI2H0jFuBRe7+s4xF9cC4MD0OuC9j/ukWOQBYl9G9Q0REWlDLs4hIeTkYOA143szmh3nfAyYBd5nZWcDrwIlh2QNEw9QtJhqq7oyCRisiUmKUPIuIlBF3/wtgrSw+Isv6DpyfaFAiImVE3TZERERERGIqq5ZnXQwgIiIiIklSy7OIiIiISExl1fIsIiJS7nTrbZF0qeVZRERERCQmJc8iIiIiIjEpeRYRERERiUnJs4iIiIhITEqeRURERERi0mgbIiIiRUyja4gUF7U8i4iIiIjEVPYtzy2/seuOgyIiIgLKESQ3ZZ88i4iIiMShZFriULcNEREREZGYlDyLiIiIiMSk5FlEREREJKZEkmcz29XMbjWze5LYv4iIiIhIGmInz2Z2m5mtNLOFLeYfY2Yvm9liM5sI4O5L3P2sfAcrIiLty1Zfm1kfM3vEzF4Jf3uH+WZmPw91+AIz+1x6kYuIFL+OtDxPBY7JnGFmXYCbgWOBIcDJZjYkb9GJiEguptKivgYmAo+5+2DgsfAcovp7cHhMAH5RoBhFREpS7KHq3H2OmdW0mL0/sNjdlwCY2QxgDPBi3iIUEZEOaaW+HgOMCNPTgNnAJWH+dHd3YK6Z9TKzndx9RYHCFSla2e7uqOHrpLPjPA8AlmU8bwSGm1lf4FpgHzO71N1/3HJDM5tA1MrBwIEDOxmGiIi0Y8eMhPgNYMcwna0eHwAoeS4A3XpbpPQkcpMUd28CzmlnnSnAFIC6ujpPIg4REdmau7uZdajeVYOHiEiks8nzcmCXjOfVYV7R0t2DRKRCvdncHcPMdgJWhvmx6nE1eIiIRDqbPD8DDDazQUSV7Vjga52OSkRE8q0eGAdMCn/vy5h/QbhmZTiwTv2dRVqnRjjpyFB1dwJPAp81s0YzO8vdNwIXAA8Bi4C73P2FZEIVEZE4stXXREnzkWb2CvCF8BzgAWAJsBj4FXBeCiGLiJSMjoy2cXIr8x8gqnxFRKQItFZfA0dkWdeB85ONSESkfOj23CIiIiIiMSUy2oaIiIhsTUPTiZQ+Jc8iIiIiOdIFhJWn4pNn3T1IREREROJSn2cRERERkZgqvuVZREQkKerjLFJ+lDyLiIiI5In6QJc/ddsQEREREYlJybOIiIiISEzqtiEiIpIH6t8sUhmUPCdEfZ5EREREQ+KWHyXPIiIiMahRRPJF51JpU59nEREREZGY1PIsIiKSA/VxFqlMSp5FREREioj6SRc3Jc8iIiIiRU79pIuHkmcREZEs1C1DCkXnWmlR8pxFeydxob7t5fKzTVI/9eSjTPStWYqJfhbdkpkdA9wIdAFucfdJKYckIlKUlDyLiFQ4M+sC3AwcCTQCz5hZvbu/mG5k+dPel3e1/EmpU4NA4Sh5FhGR/YHF7r4EwMxmAGOAkkye4yTCSpalEuTypVEJd/uUPIuIyABgWcbzRmB4SrG0S4mvSG7/B/n4YqnkGszd044BM1sFvJ7Dpv2At/IcTilSOagMmqkcIoUuh0+7e/8Cvl5emdkJwDHufnZ4fhow3N0vyFhnAjAhPP0s8HI7u9W5+DGVxcdUFltSeXyskGXRqTq7KFqecz0AM2tw97p8x1NqVA4qg2Yqh4jKocOWA7tkPK8O8zZz9ynAlLg71HvwMZXFx1QWW1J5fKyUykK35xYRkWeAwWY2yMy6AWOB+pRjEhEpSkXR8iwiIulx941mdgHwENFQdbe5+wsphyUiUpRKPXmO/RNimVM5qAyaqRwiKocOcvcHgAfyuEu9Bx9TWXxMZbEllcfHSqYsiuKCQRERERGRUqA+zyIiIiIiMaWaPJvZMWb2spktNrOJWZZ3N7OZYflTZlaTsezSMP9lMzu6vX2GC2GeCvNnhotiikKBy2Gqmb1mZvPDozbp44sroXK4zcxWmtnCFvvqY2aPmNkr4W/vRA8upgKXwVVmtjzjXBiV6MF1QL7Lwcx2MbPHzexFM3vBzC7KWL8oz4VyYGZfDeX9kZmVxFX0+dbeuVwpWquHKlFb9VGlMbMqM3vazJ4LZfGDtGOKxd1TeRBdlPIqsCvQDXgOGNJinfOAyWF6LDAzTA8J63cHBoX9dGlrn8BdwNgwPRk4N61jT7kcpgInpH3chSiHsOww4HPAwhb7+gkwMUxPBK6rwDK4Crg47eMuRDkAOwGfC+v0BP6R8T9RdOdCuTyAPYjGhJ4N1KUdTwrH3+65XCmP1uqhSny0VR9V2gMwoEeY7go8BRyQdlztPdJsed58O1h3/xBovh1spjHAtDB9D3CEmVmYP8PdP3D314DFYX9Z9xm2OTzsg7DP45M7tA4pWDkU4Fg6I4lywN3nAKuzvF7mvorlfCh0GRSrvJeDu69w92cB3P0dYBHRXfVa7qtYzoWy4O6L3L29m6mUs1KsixNRgvVQYtqpjyqKR94NT7uGR9FfjJdm8pztdrAtT57N67j7RmAd0LeNbVub3xdYG/bR2mulpZDl0OxaM1tgZjeYWfd8HEQeJFEObdnR3VeE6TeAHXMLO68KXQYAF4Rz4bYi6q6QaDmELh77ELVwQHGeC1Iecv2/lAqRpT6qOGbWxczmAyuBR9y96MtCFwxWnkuB3YH9gD7AJemGkz6Pfi8q+m+6CfgF8O9ALbAC+H+pRlMAZtYD+B3wTXd/u+XyCj4XcmZmj5rZwiyPimxhFYmrvfqoUrj7JnevJbqz6f5mNjTlkNqV5jjP7d4ONmOdRjPbFtgBaGpn22zzm4BeZrZtaKXK9lppKWQ5kNHC9oGZ3Q5cnIdjyIekyqE1b5rZTu6+wsx2IvrGm7aCloG7v9k8bWa/Av6Yc+T5lUg5mFlXog+q37r7vRnrFOO5UDLc/Qtpx1DEcqmbpAK0UR9VLHdfa2aPA8cARX1haZotz3FuB1sPjAvTJwCzQstQPTA2XHE/CBgMPN3aPsM2j4d9EPZ5X4LH1hEFKweAkBwQ+oceT/GcoEmUQ1sy91Us50NBy6D5XAi+TBmfC+F8vxVY5O4/a2NfxXIuSHnQbc9lK+3URxXFzPqbWa8wvR1wJPBSqkHFkebVisAooqtMXwUuC/OuBkaH6SrgbqKLfp4Gds3Y9rKw3cvAsW3tM8zfNexjcdhn9zSPPcVymAU8T5Qo/YZwlWsxPBIqhzuJuiRsIOpveFaY3xd4DHgFeBTok/bxp1AGvw7nwgKiD/Sd0j7+pMoBOISoO8YCYH54jCrmc6EcHkRfyhqBD4A3gYfSjimFMshaF1fao7V6qBIfbdVHlfYAhgF/D2WxELgi7ZjiPHSHQRERERGRmHTBoIiIiIhITEqeRURERERiUvIsIiIiIhJTmkPVbdavXz+vqalJOwwRkZzMmzfvLXfvn3YchaI6W0RKWWfr7KJInmtqamhoaEg7DBGRnJjZ62nHUEiqs0WklHW2zla3DRERERGRmJQ8i4iIiIjEVBTdNkQKrWbi/Vs8XzrpuJQiEREpHy3rVlD9KuVHybOIdMqGDRtobGxk/fr1aYeSuKqqKqqrq+natWvaobTKzKqAOUB3ojr+Hne/MtyufAbRHRXnAae5+4dm1h2YDuwLNAEnufvSVIIXkURUUj2dKak6W8mziHRKY2MjPXv2pKamBjNLO5zEuDtNTU00NjYyaNCgtMNpywfA4e7+rpl1Bf5iZn8CvgXc4O4zzGwycBbwi/B3jbt/xszGAtcBJ6UVvBQ3/WpXmiqlns6UZJ2t5FlEOmX9+vUVUSGbGX379mXVqlVph9Imd3fg3fC0a3g4cDjwtTB/GnAVUfI8JkwD3APcZGYW9iOSd0rAC69S6ulMSdbZSp5FpNMqpUIuleM0sy5EXTM+A9wMvAqsdfeNYZVGYECYHgAsA3D3jWa2jqhrx1sFDVrKVrZ+0B3dRgl255VK/ZVPSR2zRtuQklcz8f4tHu0tz6Uil9LWo0cPAP75z39ywgkntLre2rVr+Z//+Z9ChZUYd9/k7rVANbA/sHtn92lmE8yswcwair31XeJRvSiFZGaceuqpm59v3LiR/v3788UvfjHr+jU1Nbz1VvQd/qCDDmpz3z/60Y/yF2gM7bY86+ITEemIfH8I57PFaeedd+aee+5pdXlz8nzeeefl7TXT5O5rzexx4ECgl5ltG1qfq4HlYbXlwC5Ao5ltC+xAVHe33NcUYApAXV2dunQIkP//dymMNOrpT3ziEyxcuJD333+f7bbbjkceeYQBAwa0ux3A3/72tzaX/+hHP+J73/terH3lQ5yW5+aLT/YGaoFjzOwAootKbnD3zwBriC46gYyLT4AbwnoiJUWt1aVl+vTpDBs2jL333pvTTjuN1157jQMPPJC99tqLyy+/fPN6S5cuZejQoQC88MIL7L///tTW1jJs2DBeeeUVJk6cyKuvvkptbS3f+c530jqcTjGz/mbWK0xvBxwJLAIeB5qb3ccB94Xp+vCcsHyW+juLSBJGjRrF/fdHn6d33nknJ5988uZlTU1NHHXUUey5556cffbZZFZDzb8erlixgsMOO4za2lqGDh3KE088wcSJE3n//fepra3llFNOKchxtNvyrItPpJgoiZWWXnjhBa655hr+9re/0a9fP1avXs348eM599xzOf3007n55puzbjd58mQuuugiTjnlFD788EM2bdrEpEmTWLhwIfPnzy/sQeTXTsC00O95G+Aud/+jmb0IzDCza4C/A7eG9W8Ffm1mi4HVwNg0gpbio/pW8m3s2LFcffXVfPGLX2TBggWceeaZPPHEEwD84Ac/4JBDDuGKK67g/vvv59Zbb91q+zvuuIOjjz6ayy67jE2bNvHee+9x6KGHctNNNxW03o51wWASF5+Y2QRgAsDAgQM7dxQiGVThV5ZZs2bx1a9+lX79+gHQp08f/vrXv/K73/0OgNNOO41LLrlkq+0OPPBArr32WhobG/mP//gPBg8eXNC4k+LuC4B9ssxfQtT/ueX89cBXCxCaFDnVnZK0YcOGsXTpUu68805GjRq1xbI5c+Zw7733AnDcccfRu3fvrbbfb7/9OPPMM9mwYQPHH388tbW1hQh7K7EuGEzi4hN3n+Lude5e179//87uTkRkC+1dZf21r32N+vp6tttuO0aNGsWsWbMKFJmISOUaPXo0F1988RZdNuI67LDDmDNnDgMGDGD8+PFMnz49gQjb16HRNtx9LVG/uc0Xn4RF2S4+oa2LT0RE8uHwww/n7rvvpqkpqmZWr17NwQcfzIwZMwD47W9/m3W7JUuWsOuuu3LhhRcyZswYFixYQM+ePXnnnXcKFruISKU588wzufLKK9lrr722mH/YYYdxxx13APCnP/2JNWvWbLXt66+/zo477sjXv/51zj77bJ599lkAunbtyoYNG5IPPmg3edbFJyJSzPbcc08uu+wyPv/5z7P33nvzrW99ixtvvJGbb76Zvfbai+XLl2fd7q677mLo0KHU1taycOFCTj/9dPr27cvBBx/M0KFDS/aCQRGRYlZdXc2FF1641fwrr7ySOXPmsOeee3Lvvfdm7dI7e/Zs9t57b/bZZx9mzpzJRRddBMCECRMYNmxYwS4YtPbyWjMbRnRBYObFJ1eb2a5EQ9X1Ibr45FR3/yAMbfdroj53q4Gxoa9dq+rq6ryhoaHTByPlL6k+eS2H2cn2OhqkP7tFixaxxx57pB1GwWQ7XjOb5+51KYVUcKqzy0Ox9HHOVrfqJin5VWn1dKYk6uw4o23o4hORLFS5i4h0XrEk8SJx6fbcUtRUqYqIiEgxUfIsEpMSeREREenQaBsiItlUyjXBlXKcIlJ+KrH+SuqY1fIsRSWt1l21KueuqqqKpqYm+vbt2+7YyqXM3WlqaqKqqirtUEREOqRS6ulMSdbZSp5FpFOqq6tpbGxk1apVaYeSuKqqKqqrq9MOQ0SkQyqpns6UVJ2t5FlSo9be8tC1a1cGDRqUdhgiUsY0ulHnqJ7OLyXPIiIiUlI0Fr+kSRcMioiIiIjEpJZnKRh10xAREZFSp+RZRESkzKnxQiR/1G1DRERERCQmJc8iIiIiIjEpeRYRERERiUnJs4iIiIhITLpgUCRPNO6oFAMz2wWYDuwIODDF3W80sz7ATKAGWAqc6O5rLLpX743AKOA9YLy7P5tG7CIipUDJs4hIedkIfNvdnzWznsA8M3sEGA885u6TzGwiMBG4BDgWGBwew4FfhL9SwjS6hkhylDyLiJQRd18BrAjT75jZImAAMAYYEVabBswmSp7HANPd3YG5ZtbLzHYK+xEpWbqltyRFybOISJkysxpgH+ApYMeMhPgNom4dECXWyzI2awzzlDxLSVFruxRKu8mz+s9JHOrvK1JczKwH8Dvgm+7+dlQ1R9zdzcw7uL8JwASAgQMH5jNUyQMljiKFE6flWf3nRERKiJl1JUqcf+vu94bZbzZ3xzCznYCVYf5yYJeMzavDvC24+xRgCkBdXV2HEm+RYqBGHsmXdpNn9Z+TXKklRH3upPDCr3+3Aovc/WcZi+qBccCk8Pe+jPkXmNkMooaOdaqvRURa16E+z/nsP6efAEVEEnEwcBrwvJnND/O+R5Q032VmZwGvAyeGZQ8QdbNbTNTV7oyCRisiUmJiJ8/57j+nnwBFRPLP3f8CWCuLj8iyvgPnJxqUiEgZiXWHwbb6z4XlHe4/JyIiIiJSauKMtqH+cyIiIkVC15OIpCtOtw31nxPJE11AKCIiUtrijLah/nMiIiIiIsTs8ywiIiIiIkqeRURERERi69A4zyIiIlJYukBQpLgoeZZYdKFbMnS7WBERkdKibhsiIiIiIjEpeRYRERERiUnJs4iIiIhITOrzLDnRBSwiIiJSiZQ8i4iIFBE1TogUNyXPIkVGI5uIiBSG6lvJhfo8i4iIiIjEpORZRERERCQmdduQrNTnTkRERGRrSp5FREREUB9oiUfdNkREREREYlLyLCIiIiISk7ptiIiUGTO7DfgisNLdh4Z5fYCZQA2wFDjR3deYmQE3AqOA94Dx7v5sGnFXAnULECl9Sp5FRMrPVOAmYHrGvInAY+4+ycwmhueXAMcCg8NjOPCL8FcKQBdnF7ds74++8Ei7ybNaMERESou7zzGzmhazxwAjwvQ0YDZR8jwGmO7uDsw1s15mtpO7ryhQuCIlRb8eSJyW56moBaOsqeVDpCLsmJEQvwHsGKYHAMsy1msM85Q8i4hk0W7yrBYMEZHy4u5uZt6RbcxsAjABYODAgYnEJVKK1BJdeXIdbaOjLRhbMbMJZtZgZg2rVq3KMQwREYnpTTPbCSD8XRnmLwd2yVivOszbgrtPcfc6d6/r379/4sGKiBSrTl8wmEsLRthuCjAFoK6ursPbi1SKON1q1NIhMdQD44BJ4e99GfMvMLMZRN3s1unXQhGR1uWaPL/Z3B0jlxYMERFJjpndSdS1rp+ZNQJXEiXNd5nZWcDrwIlh9QeILvJeTHSh9xkFD1hEpITkmjyrBUNEpEi5+8mtLDoiy7oOnJ9sRJVLF2SLlJ84Q9WpBUNEREREhHijbagFo8yoJUREREQkN7rDoEgZ0FBJIiIihaHkuQKopVlEREQkP5Q8lxklyiIiIiLJyfUmKSIiIiIiFUctzyJlSH2gRZKn/zOJI9svwjpXSpuSZxERkTxQtzkBnQeVQMmzSAVQy4eIiEh+KHkWERHJQt0yRCQbXTAoIiIiIhKTWp5FKpRa1URERDpOyXOJ04UJIiKFofpW8kWNF6VN3TZERERERGJSy3OJUcuHiIhIeVFLdGlRy7OIiIiISExqeS5iamWWQlLLh1QSne9SzDQ2f3FT8iwiWcWpvJWASKlQY4SI5IuSZxGJTQmIlIJczlOd2yISl5LnFKmyFhERkTj0S1/xSCx5NrNjgBuBLsAt7j4pqdcSkeKgyr10qc4WKS2qb9OTSPJsZl2Am4EjgUbgGTOrd/cXk3i9YtBeK7JOaqlEuuilNJR6na1f8UTi/R+o/s2PpFqe9wcWu/sSADObAYwBSqIibo/604nkT5zWE7WwJC61OlvvrUjh5JJg6390a0klzwOAZRnPG4HhmSuY2QRgQnj6rpm93InX6we81Ynt803xtE3xtK2s47HrOrcc6GfXFV35fDrtIDqpaOrsGO9/PhXT/5piya5YYimWOCDhWDpYRxdLuXQ0jk7V2aldMOjuU4Ap+diXmTW4e10+9pUPiqdtiqdtiqdtRRpPTdpxJK0c6+xiiQMUS2uKJZZiiQMUSzHEkdQdBpcDu2Q8rw7zRESk+KjOFhGJKank+RlgsJkNMrNuwFigPqHXEhGRzlGdLSISUyLdNtx9o5ldADxENOzRbe7+QhKvFeTlp8Q8UjxtUzxtUzxtUzx5VsF1drHEAYqlNcUSS7HEAYolm4LGYe5eyNcTERERESlZSXXbEBEREREpO0qeRURERERiKpnk2cz6mNkjZvZK+Ns7yzq1Zvakmb1gZgvM7KSMZVPN7DUzmx8etSnHM8jMnjKzxWY2M1ykk2g8Yb0HzWytmf2xxfyCl0878aRVPuPCOq+Y2biM+bPN7OWM8vlUjnEcE/az2MwmZlnePRzv4nD8NRnLLg3zXzazo3N5/XzEYmY1ZvZ+RllM7mwsMeM5zMyeNbONZnZCi2VZ37cU49mUUT4VeeGdmf3UzF4Kdd/vzaxXxrJ2z+V81QFm9tVQB39kZnUZ87uZ2e1m9ryZPWdmI1rZvtbM5ob3ssHM9s8ljnzEEtb9RijXF8zsJ2nGEtb/tpm5mfVLI462zrMUYon1OdPJWLqa2bQQyyIzu7SV7Y8I9dN8M/uLmX0mpTjMzK41s3+E9S7MJY58xJKx/s/N7N1c4wDA3UviAfwEmBimJwLXZVlnN2BwmN4ZWAH0Cs+nAicUUTx3AWPD9GTg3KTjCcuOAL4E/LHF/IKXTzvxFLx8gD7AkvC3d5juHZbNBuo6GUMX4FVgV6Ab8BwwpMU65wGTw/RYYGaYHhLW7w4MCvvpklIsNcDCfJ0rHYinBhgGTM88V9t639KIJyx7N5/lU4oP4Chg2zB9XfP/XNxzOV91ALAH8NmW/8PA+cDtYfpTwDxgmyzbPwwcG6ZHAbM7USadjWUk8CjQvXndtGIJy3chusj0daBfSmWS9TxLKZZYn3udjOVrwIwwvT2wFKjJsv0/gD3C9HnA1JTiOIOojtwmwXM2VixheR3wazpZR5dMyzPRrWKnhelpwPEtV3D3f7j7K2H6n8BKoH+xxWNmBhwO3NPW9vmOJ8TxGPBOJ18r0XhSLJ+jgUfcfbW7rwEeAY7p5Otm2nwLZHf/EGi+BXJrcd4DHBHKYwxR5fCBu78GLA77SyOWJLQbj7svdfcFwEcttk3ifetMPAK4+8PuvjE8nUs0djTEOJfzWQe4+yJ3z3Y3xCHArLDOSmAt0QfrVrsAPhmmdwD+mUsceYrlXGCSu3+QsW5asQDcAHyXqIxSiaON86zgsRDzc6+TsTjwCTPbFtgO+BB4u5X1On3e5iGOc4Gr3f2jsL8kztlYsZhZF+CnROdsp5RS8ryju68I028AO7a1cvhprRtRq0aza8NPOzeYWfcU4+kLrM34h28kuj1uweJpRWrl00Ja5ZPtFsWZr3t7+Ans+zkmke3tf4t1wvGvIyqPONsWKhaAQWb2dzP7s5kd2ok4OhJPEtsmtc8qi37in2tmx3cylnJwJvCnMB2nbJOoA1p6DhhtZtua2SBgX7a8UUyzbwI/NbNlwPVAmz8JJxzLbsChFnVn+bOZ7ZdWLGY2Blju7s8lEEPsOFrIPM/SiCUfn8PtuQf4F9Ev2f8HXO/uq7OsdzbwgJk1AqcBk1KK49+Bk0J9+CczG5znODoSywVAfcZ7lLPUbs+djZk9CvxblkWXZT5xdzezVr/pmtlORM3y45q/7RBVeG8QJbBTgEuAq9OIJ9fGu3zF04rUyidfEo7nFHdfbmY9gd8RVUbTc4u05K0ABrp7k5ntC/zBzPZ092ytDpXq0+F82RWYZWbPu/ur7W5VYtr6n3P3+8I6lwEbgd+mGUcWtxH9DNxA1O3gb8CmLOudC/yXu//OzE4EbgW+kFIs2xJ1TzoA2A+4y8x29fB7dKFiMbPtge8RdZloV8Jl0vwasc6zQsQC8T5ncoxl//DaOxN1U3vCzB519yUt1vsvYJS7P2Vm3wF+RpRQFzqO7sB6d68zs/8gKstWG1ySisXMdga+Coxo7bU7oqiSZ3dvq0J608x2cvcVIRnN2vRvZp8E7icq6LkZ+27+pvGBmd0OXJxiPE1ALzPbNrSsxLoVbj7iaWPfqZRPK9Iqn+Vs+Y9VTdS3CndfHv6+Y2Z3EP2zdjR5jnML5OZ1GsNPUDsQlUe+b5+ccyzhg7r5Z+N5ZvYqUYtYQ8LxtLXtiBbbzu5ELJ2NJ/N8WWJms4F92PJXsLLQ1v8cgJmNB74IHJGR4MUp2w7VAe3F0co2G4kSjOZY/0bUT7SlccBFYfpu4JZ29ptkLI3AvaEsnzazj4B+wKoCx/LvRP3VnwuNQdXAs2a2v7u/UcA4mpeNZ+vzrLX9JhlLhz73comFqH/vg+6+AVhpZn8l6kKSmSj2B/Z296fCrJnAg4WOI2gE7g3Tvwdub2unCcayD/AZYHE4Z7c3s8XuntOFlKXUbaOeqBIj/N3qG4hFV2P/Hpju7ve0WLZT+GtE/ZAWphVP+Od+HDihre3zHU9b0iif1qRYPg8BR5lZb4uukj4KeCj8VNcPoqt6iSrpXMonzi2QM+M8AZgVyqMeGGvRCBiDgMHA0znE0OlYzKy/RX3HCC2rg9m6wkwintZkfd/SiifE0T1M9wMOBl7sZDwlx8yOIepbONrd38tY1O65nFAd0DK+7c3sE2H6SGCju2d7n/4JfD5MHw68ks84OhjLH4guGsTMdiP6pfCtQsfi7s+7+6fcvcbda4gSpM9lS5yTjCMsa+08y5sOvD+d+hyO6f+IzkNCTAcAL7VYZw2wQzhHAI4EFqUQB2Scs0T/R1m/ACUdi7vf7+7/lnHOvpdr4ty8w5J4EPWBe4yo4noU6OMfXzl5S5g+FdgAzM941IZls4DniZKe3wA9Uo5nV6IPjMVErRndk44nPH+CqJXifaIK7+i0yqedeNIqnzPDay4GzgjzPkF0dfUC4AXgRnIc6YLoav1/ELVCXhbmXU1U8QNUheNdHI5/14xtLwvbvUy4+r+TZZJTLMBXQjnMB54FvtTZWGLGs184R/5F1DL5QlvvW1rxAAeF/6Xnwt+z8hFPqT3Ce7GMj+u+ye2dy8ADwM5hOi91APDl8D59ALwJPBTm14TXX0RUJ3w6Y5tbCFfzA4eE///ngKeAfTtRJp2NpRtR/bww/O8dnlYsLfa1lNxH2+hsmbR6nqUQS9bPmTzH0iP8P7xA9KX8O638/3yZj+uh2WR8lhQ4jl5Ev74/DzxJ1CKeSpm02FenRtvQ7blFREqcmd1G9IvISncfmmW5EX3pGwW8B4x392fDsnHA5WHVa9x9WsvtRUTkY6XUbUNERLKbStvD8x1L1D1iMDAB+AVEN3UArgSGE/Xjv9I6cWMHEZFKoORZRKTEufscINvQTM3GEF174R5duNwrXOeQ9NjmIiJlpyhG2+jXr5/X1NSkHYaISE7mzZv3lrsndUOmfGhtbOXY41mb2QSiVms+8YlP7Lv77rsnE6mISMI6W2cXRfJcU1NDQ0NnRrkSEUmPmb2edgxJc/cpRGPAU1dX56qzRaRUdbbOVrcNEZHy19rYyvkeP1xEpOwpeRYRKX/1wOkWOQBY59GNkZIYI1tEpKwVRbcNERHJnZndSXSXxX5m1kg0gkZXAHefTDTW6SiiMXHfA84Iy1ab2Q+JbgoDcLW7t3XhoYhIxSva5HnDhg00Njayfv36tEMpGlVVVVRXV9O1a9e0QxFJTM3E+7eat3TScSlEUjrc/eR2ljtwfivLbgNuSyIuEZFyVLTJc2NjIz179qSmpoZwH/KK5u40NTXR2NjIoEGD0g5HREREpCIVbZ/n9evX07dvXyXOgZnRt29ftcSLiIiIpKhok2dAiXMLKg8RERGRdBVttw0RKU8t+zTH6c+cyzYiIiJJKJnkOdtFRJ2Rrw/fq666ih49enDxxRe3ud748eP54he/yAknnMDZZ5/Nt771LYYMGZJ13alTp3LUUUex88475yVGEREREcmPkkmey8ktt9zS5vKpU6cydOhQJc9SEfL9xVhERCRJRd3nOW1Lly5l9913Z/z48ey2226ccsopPProoxx88MEMHjyYp59+GoDnnnuOAw88kMGDB/OrX/0KiEbHuOCCC/jsZz/LF77wBVauXLl5vyNGjKChoYFNmzYxfvx4hg4dyl577cUNN9zAPffcQ0NDA6eccgq1tbW8//77qRy7iIiIiGxNLc/tWLx4MXfffTe33XYb++23H3fccQd/+ctfqK+v50c/+hG1tbUsWLCAuXPn8q9//Yt99tmH4447jrlz5/Lyyy/z4osv8uabbzJkyBDOPPPMLfY9f/58li9fzsKFCwFYu3YtvXr14qabbuL666+nrq4ujUMWERERkVao5bkdgwYNYq+99mKbbbZhzz335IgjjsDM2GuvvVi6dCkAY8aMYbvttqNfv36MHDmSp59+mjlz5nDyySfTpUsXdt55Zw4//PCt9r3rrruyZMkSvvGNb/Dggw/yyU9+ssBHJyIiIiIdoeS5Hd27d988vc0222x+vs0227Bx40Zg6yHk4g4p17t3b5577jlGjBjB5MmTOfvss/MUtYiIiIgkIXbybGZdzOzvZvbH8HyQmT1lZovNbKaZdQvzu4fni8PymoRiLxr33Xcf69evp6mpidmzZ7Pffvtx2GGHMXPmTDZt2sSKFSt4/PHHt9rurbfe4qOPPuIrX/kK11xzDc8++ywAPXv25J133in0YYiIiIhIOzrS5/kiYBHQ3LfgOuAGd59hZpOBs4BfhL9r3P0zZjY2rHdSZwMt5nFdhw0bxsiRI3nrrbf4/ve/z84778yXv/xlZs2axZAhQxg4cCAHHnjgVtstX76cM844g48++giAH//4x0A0rN0555zDdtttx5NPPsl2221X0OMRERERkezM3dtfyawamAZcC3wL+BKwCvg3d99oZgcCV7n70Wb2UJh+0sy2Bd4A+nsbL1RXV+cNDQ1bzFu0aBF77LFHrsdVtlQuUuryMTRdsX2ZNrN57l4xV/hmq7NFREpFZ+vsuN02/j/gu8BH4XlfYK27bwzPG4EBYXoAsAwgLF8X1t+CmU0wswYza1i1alVu0YuIiIiIFFC73TbM7IvASnefZ2Yj8vXC7j4FmAJRK0a+9isi6dENT9JjZscANwJdgFvcfVKL5TcAI8PT7YFPuXuvsGwT8HxY9n/uProgQYuIlKA4fZ4PBkab2SigiqjP841ALzPbNrQuVwPLw/rLgV2AxtBtYwegKZfg3D32yBWVIE4XGxGpPGbWBbgZOJLol8BnzKze3V9sXsfd/ytj/W8A+2Ts4n13ry1QuCIiJa3dbhvufqm7V7t7DTAWmOXupwCPAyeE1cYB94Xp+vCcsHxWW/2dW1NVVUVTU5MSxsDdaWpqoqqqKu1QRKT47A8sdvcl7v4hMAMY08b6JwN3FiQyEZEy05k7DF4CzDCza4C/A7eG+bcCvzazxcBqooS7w6qrq2lsbET9oT9WVVVFdXV12mGISPHZfK1J0AgMz7aimX0aGATMyphdZWYNwEZgkrv/IaE4RURKXoeSZ3efDcwO00uIWjtarrMe+GpnA+vatSuDBg3q7G5ERGRLY4F73H1TxrxPu/tyM9sVmGVmz7v7q5kbmdkEYALAwIEDCxetiEiR6UzLs4hUOF0gWDSarzVplnkdSktjgfMzZ7j78vB3iZnNJuoP/WqLdXSRt4gIuj23iEg5eAYYHO782o0oQa5vuZKZ7Q70Bp7MmNfbzLqH6X5EF4m/2HJbERGJqOVZRKTEhZtVXQA8RDRU3W3u/oKZXQ00uHtzIj0WmNHiIu49gF+a2UdEDSqTMkfpEBGRLSl5FhEpA+7+APBAi3lXtHh+VZbt/gbslWhwIiJlRN02RERERERiUsuziJSclhcqLp10XEqRiIhIpVHLs4iIiIhITEqeRURERERiUvIsIiIiIhKT+jyLSGy6KYqIiFQ6tTyLiIiIiMSk5FlEREREJCYlzyIiIiIiMSl5FhERERGJScmziIiIiEhMSp5FRERERGJS8iwiIiIiEpOSZxERERGRmJQ8i4iIiIjEpDsMikhWuptgaTGzY4AbgS7ALe4+qcXy8cBPgeVh1k3ufktYNg64PMy/xt2nFSRoEZESpORZRAAly6XMzLoANwNHAo3AM2ZW7+4vtlh1prtf0GLbPsCVQB3gwLyw7ZoChC4iUnLUbUNEpPTtDyx29yXu/iEwAxgTc9ujgUfcfXVImB8BjkkoThGRkqfkWUSk9A0AlmU8bwzzWvqKmS0ws3vMbJcObisiIqjbhoiUgWxdTpZOOi6FSIra/wJ3uvsHZvafwDTg8Lgbm9kEYALAwIEDk4lQRKQEKHkWqVDq41xWlgO7ZDyv5uMLAwFw96aMp7cAP8nYdkSLbWe3fAF3nwJMAairq/POBiwiUqrUbUNEpPQ9Aww2s0Fm1g0YC9RnrmBmO2U8HQ0sCtMPAUeZWW8z6w0cFeaJiEgWankWESlx7r7RzC4gSnq7ALe5+wtmdjXQ4O71wIVmNhrYCKwGxodtV5vZD4kScICr3X11wQ9CRKREKHkWESkD7v4A8ECLeVdkTF8KXNrKtrcBtyUaoIhImWg3eQ5XZE8HdiQaA3SKu98YxgadCdQAS4ET3X2NmRnRQP2jgPeA8e7+bDLhi0gc6t8sIiKSH3H6PG8Evu3uQ4ADgPPNbAgwEXjM3QcDj4XnAMcCg8NjAvCLvEctIiIiIpKCdpNnd1/R3HLs7u8QXWQygGgA/uZbuE4Djg/TY4DpHpkL9GpxoYqIiIiISEnq0GgbZlYD7AM8Bezo7ivCojeIunWABtwXERERkTIVO3k2sx7A74Bvuvvbmcvc3Yn6Q8dmZhPMrMHMGlatWtWRTUVEREREUhEreTazrkSJ82/d/d4w+83m7hjh78owv93B+iEacN/d69y9rn///rnGLyIiIiJSMO0mz2H0jFuBRe7+s4xF9cC4MD0OuC9j/ukWOQBYl9G9Q0RERESkZMUZ5/lg4DTgeTObH+Z9D5gE3GVmZwGvAyeGZQ8QDVO3mGioujPyGbBIpYsz7NzSSccVIBIREZHK027y7O5/AayVxUdkWd+B8zsZl4iIiIhI0enQaBsiIiIiIpVMybOIiIiISExx+jyLSInR7bhFRESSoZZnEREREZGYlDyLiIiIiMSk5FlEpAyY2TFm9rKZLTaziVmWf8vMXjSzBWb2mJl9OmPZJjObHx71hY1cRKS0qM+zSJFT/2Vpj5l1AW4GjgQagWfMrN7dX8xY7e9Anbu/Z2bnAj8BTgrL3nf32kLGLCJSqtTyLCJS+vYHFrv7Enf/EJgBjMlcwd0fd/f3wtO5QHWBYxQRKQtqeRaRstSyxb7M77o4AFiW8bwRGN7G+mcBf8p4XmVmDcBGYJK7/yHvEYqIlAklzyIpytYlo8yTPEmZmZ0K1AGfz5j9aXdfbma7ArPM7Hl3f7XFdhOACQADBw4sWLwiIsVG3TZERErfcmCXjOfVYd4WzOwLwGXAaHf/oHm+uy8Pf5cAs4F9Wm7r7lPcvc7d6/r375/f6EVESoiSZxGR0vcMMNjMBplZN2AssMWoGWa2D/BLosR5Zcb83mbWPUz3Aw4GMi80FBGRDOq2IVJkNLqGdJS7bzSzC4CHgC7Abe7+gpldDTS4ez3wU6AHcLeZAfyfu48G9gB+aWYfETWoTGoxSoeIiGRQ8iwiUgbc/QHggRbzrsiY/kIr2/0N2CvZ6EREyoe6bYiIiIiIxKTkWUREREQkJnXbEEmQ+i+LiIiUF7U8i4iIiIjEpORZRERERCQmddsQyRN10RARESl/ankWEREREYlJybOIiIiISExKnkVEREREYlKfZ5EssvVfXjrpuBQiERERkWKi5FkkJl0QWNpavn/6MiQiIrlQtw0RERERkZjU8iwVSa3IIiIikgu1PIuIiIiIxKSWZyk7alUWERGRpCSWPJvZMcCNQBfgFneflNRrSeVQYiySXXt1rpl1B6YD+wJNwEnuvjQsuxQ4C9gEXOjuDxUwdBGRkpJI8mxmXYCbgSOBRuAZM6t39xeTeD0pTu2NbqBEWCQ/Yta5ZwFr3P0zZjYWuA44ycyGAGOBPYGdgUfNbDd331TYoxARKQ1JtTzvDyx29yUAZjYDGAMoeS5R+RjmS8myFJMyG8s7Tp07BrgqTN8D3GRmFubPcPcPgNfMbHHY35MFil1EpKQklTwPAJZlPG8Ehif0WgVXqA/dpF4nH0msEmGRohKnzt28jrtvNLN1QN8wf26LbQckF6qISGlL7YJBM5sATAhP3zWzl9OKpRX9gLfirmzXJRhJcq/ToWMsQeV+fKBjzKtO/H99Oo9hFKUWdfYHZrYwzXhSUAn/ay3pmCtDJR7zZzuzcVLJ83Jgl4zn1WHeZu4+BZiS0Ot3mpk1uHtd2nEkqdyPsdyPD3SMslm7dW7GOo1mti2wA9GFg3G23aLOrsT3RMdcGXTMlcHMGjqzfVLjPD8DDDazQWbWjehilPqEXktEpNLFqXPrgXFh+gRglrt7mD/WzLqb2SBgMPB0geIWESk5ibQ8h/50FwAPEQ2bdJu7v5DEa4mIVLrW6lwzuxpocPd64Fbg1+GCwNVECTZhvbuILi7cCJyvkTZERFqXWJ9nd38AeCCp/RdA0XYpyaNyP8ZyPz7QMUqQrc519ysyptcDX21l22uBazvwcpX4nuiYK4OOuTJ06pgt+tVORERERETak1SfZxERERGRsqPkGTCz28xsZebQS2ZWa2ZzzWy+mTWY2f5pxtgZZraLmT1uZi+a2QtmdlGY38fMHjGzV8Lf3mnHmqs2jvGnZvaSmS0ws9+bWa+UQ81Ja8eXsfzbZuZm1i+tGDurrWM0s2+E9/EFM/tJmnFWEjM7xsxeNrPFZjYxy/LuZjYzLH/KzGpSCDOvYhzzt8I5usDMHjOzkh+msL1jzljvK6GeKfmRGeIcs5mdmFEf3VHoGPMtxrk9MNTBfw/n96g04syXbLldi+VmZj8P5bHAzD4Xe+fuXvEP4DDgc8DCjHkPA8eG6VHA7LTj7MTx7QR8Lkz3BP4BDAF+AkwM8ycC16UdawLHeBSwbZh/XakeY2vHF57vQnSh2OtAv7RjTeA9HAk8CnQPyz6VdqyV8CC68PBVYFegG/Bc8zmXsc55wOQwPRaYmXbcBTjmkcD2YfrcSjjmsF5PYA7RDXXq0o67AO/zYODvQO/wvKTrnZjHPAU4N0wPAZamHXcnj3mr3K7F8lHAnwADDgCeirtvtTwD7j6H6OrzLWYDnwzTOwD/LGhQeeTuK9z92TD9DrCI6A5iY4BpYbVpwPGpBJgHrR2juz/s7hvDanOJxrAtOW28hwA3AN8lOmdLVhvHeC4wyaPbR+PuK9OLsqJsvuW3u38INN/yO1NmHXIPcISZWQFjzLd2j9ndH3f398LTkq1TMsR5nwF+SNQAsb6QwSUkzjF/HbjZ3ddAWdQ7cY65bPIeaDW3yzQGmO6RuUAvM9spzr6VPLfum8BPzWwZcD1wabrh5Ef4WXUf4ClgR3dfERa9AeyYVlz51OIYM51J9C2zpGUen5mNAZa7+3PpRpVfLd7D3YBDQ7eAP5vZfqkGVzmy3fK75W27t7jlN9B8y+9SFeeYM51F6dcp7R5z+Dl7F3e/v5CBJSjO+7wbsJuZ/TV04TymYNElI84xXwWcamaNRCP3fKMwoaWmo//vmyl5bt25wH+5+y7AfxGNkVrSzKwH8Dvgm+7+duYyj37DKOmWS2j9GM3sMqIxbH+bVmz5kHl8RMfzPeCKtrYpNVnew22BPkQ/q30HuKvEWzelDJjZqUAd8NO0Y0mSmW0D/Az4dtqxFNi2RF03RgAnA78q1WtmOuBkYKq7VxN1afh1eP+lBRVK68YB94bpu4l+8ihZZtaVKCH5rbs3H9ebzT9RhL8l/bNUK8eImY0HvgicEr4klKQsx/fvwCDgOTNbSvTz8bNm9m/pRdk5rbyHjcC94ae1p4GPgJK9MLKEdOSW39iWt/wuVbFuVW5mXwAuA0Y3dycqYe0dc09gKDA71DMHAPUlftFgnPe5Eah39w3u/hrRNRiDCxRfEuIc81nAXQDu/iRQRXnXtbH+37NR8ty6fwKfD9OHA6+kGEunhFa6W4FF7v6zjEWZt+sdB9xX6NjypbVjDD+1fZfoQ+691rYvdtmOz92fd/dPuXuNu9cQVfafc/c3Ugw1Z22cp38gukgLM9uN6GKXtwoeYOXpzC2/S1W7x2xm+wC/JKpTSrrBIWjzmN19nbv3y6hn5hIde0M64eZFnHP7D0Stzlg0itFuwJICxphvcY75/4AjAMxsD6LkeVVBoyyseuD0MOrGAcC6jK6sbUvrKshiegB3AiuADUQJyFnAIcA8oitSnwL2TTvOThzfIURdMhYA88NjFFHfxMeIvhg8CvRJO9YEjnExUZ+m5nmT0441n8fXYp2llPZoG629h92A3wALgWeBw9OOtVIeofz/QXSV/mVh3tVEyRNEH653h/+zp4Fd0465AMf8KPBmxjlan3bMSR9zi3VnU+KjbcR8n42ou8qLwPPA2LRjLsAxDwH+GvKe+cBRacfcyePNltudA5yT8R7fHMrj+Y6c17rDoIiIiIhITOq2ISIiIiISk5JnEREREZGYlDyLiIiIiMSk5FlEREREJCYlzyIiIiIiMSl5FhERERGJScmziIiIiEhMSp5FRERERGL6/wEYTFdvmJanjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(4, 2, figsize=(12, 8))\n",
    "for i, name in enumerate(names):\n",
    "    row, col = divmod(i, 2)\n",
    "    axes[row, col].hist(params[:, i], bins=64, label=name)\n",
    "    axes[row, col].legend()\n",
    "    if name == \"x0dist\":\n",
    "        # axes[row, col].set_xscale('log')\n",
    "        axes[row, col].set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "7b528107-64fa-475a-8d40-0324e9a58f62",
   "metadata": {},
   "outputs": [],
   "source": [
    "#保存模拟数据的文件名格式\n",
    "def random_name():\n",
    "    return dt_.datetime.now().strftime(\"%Y%m%d_%H%M%S_\") + \\\n",
    "            str(rnd.randint(0, 1000000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75bcfa60-ab01-41a7-862b-b30eee69fd62",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "模拟随机光变:   0% 1/4675 [00:00<13:32,  5.75it/s]/usr/local/lib/python3.8/dist-packages/sncosmo/fitting.py:187: RuntimeWarning: Dropping following bands from data: csst_u, csst_NUV(out of model wavelength range)\n",
      "  warnings.warn(\"Dropping following bands from data: \" +\n",
      "模拟随机光变:   0% 2/4675 [00:22<17:18:54, 13.34s/it]/usr/local/lib/python3.8/dist-packages/sncosmo/fitting.py:187: RuntimeWarning: Dropping following bands from data: csst_u, csst_g(out of model wavelength range)\n",
      "  warnings.warn(\"Dropping following bands from data: \" +\n",
      "模拟随机光变:   0% 6/4675 [01:20<23:55:13, 18.44s/it]/usr/local/lib/python3.8/dist-packages/sncosmo/fitting.py:187: RuntimeWarning: Dropping following bands from data: csst_NUV(out of model wavelength range)\n",
      "  warnings.warn(\"Dropping following bands from data: \" +\n",
      "模拟随机光变:   0% 14/4675 [02:57<14:04:00, 10.86s/it]/usr/local/lib/python3.8/dist-packages/sncosmo/fitting.py:187: RuntimeWarning: Dropping following bands from data: csst_NUV, csst_u(out of model wavelength range)\n",
      "  warnings.warn(\"Dropping following bands from data: \" +\n",
      "模拟随机光变:   1% 26/4675 [05:46<16:47:13, 13.00s/it]/usr/local/lib/python3.8/dist-packages/sncosmo/fitting.py:187: RuntimeWarning: Dropping following bands from data: csst_g, csst_u(out of model wavelength range)\n",
      "  warnings.warn(\"Dropping following bands from data: \" +\n",
      "模拟随机光变:   1% 34/4675 [08:27<30:10:42, 23.41s/it]/usr/local/lib/python3.8/dist-packages/sncosmo/fitting.py:187: RuntimeWarning: Dropping following bands from data: csst_u(out of model wavelength range)\n",
      "  warnings.warn(\"Dropping following bands from data: \" +\n",
      "模拟随机光变:   1% 67/4675 [17:20<18:37:01, 14.54s/it]"
     ]
    }
   ],
   "source": [
    "#光变曲线模拟和拟合，并把success的SN记录下来\n",
    "good = 0\n",
    "bad = 0\n",
    "bad_input = 0\n",
    "\n",
    "for sn_z, sn_t0, sn_x0, sn_x1, sn_c in tqdm.tqdm(params[:,0:5], desc=\"模拟随机光变\"):\n",
    "    lcs = lc_simu(sn_z=sn_z, \n",
    "                  sn_t0=sn_t0, \n",
    "                  sn_x0=sn_x0, \n",
    "                  sn_x1=sn_x1, \n",
    "                  sn_c=sn_c, \n",
    "                  csst_mjd_start=56150, \n",
    "                  csst_mjd_duration=365.25*2, \n",
    "                  csst_mjd_cadence= 10, \n",
    "                  csst_exposure_time=150)\n",
    "    # lc_simu之后, 注册表中的csst_band已被裁减\n",
    "\n",
    "    # 恢复注册表中的csst_band为原始值\n",
    "    reset_csst_filter()\n",
    "    try:\n",
    "        # 判断至少某个band中SNR>3的点不少于1个\n",
    "        snr_mask = lcs[0]['flux'] / lcs[0]['fluxerr'] >= 5\n",
    "        is_bad_input = (pd.Series(lcs[0]['band'][snr_mask]).value_counts() >= 3).sum() == 0\n",
    "        if is_bad_input:\n",
    "            bad_input += 1\n",
    "            continue\n",
    "        \n",
    "        \n",
    "        result, fitted_model = sncosmo.fit_lc(\n",
    "            lcs[0], \n",
    "            model, \n",
    "            ['z', 't0', 'x0', 'x1', 'c'],  # parameters of model to vary\n",
    "            bounds={'z':(max(sn_z-0.05,0),sn_z+0.05)}\n",
    "        #         print(\"good\")\n",
    "        )\n",
    "\n",
    "        #bounds={'z':(0.0,1.3),'t0':(sn_t0-10, sn_t0+10),'x1':(-3,3),'c':(-0.3,0.3)}\n",
    "        #         print(\"good\")\n",
    "\n",
    "#         # ============================================\n",
    "#         #      老版本代码块: 仅记录拟合成功的SN\n",
    "#         #  \n",
    "#         # 添加、取消注释的快捷键: \n",
    "#         #      选中 === 到 --- 的代码块, \n",
    "#         #      按下　command + / \n",
    "#         #\n",
    "        \n",
    "#         if result.success:\n",
    "#             good +=1   \n",
    "#             # 仅当sn fit成功时，记录光变曲线\n",
    "#             with open(\"fitted_\" + random_name(), \"wb\") as FILE:\n",
    "#                 pickle.dump({\"result\":result, \n",
    "#                              \"fitted_model\":fitted_model, \n",
    "#                              \"lcs\":lcs, \n",
    "#                              \"true_values\":[sn_z, sn_t0, sn_x0, sn_x1, sn_c]}, file=FILE)\n",
    "#         else:\n",
    "#             bad += 1\n",
    "#         #\n",
    "#         # ---------------------------------------------\n",
    "\n",
    "        # ============================================\n",
    "        #      新版本代码块: 记录所有参与拟合的SN\n",
    "        #  \n",
    "        # 添加、取消注释的快捷键: \n",
    "        #      选中 === 到 --- 的代码块, \n",
    "        #      按下　command + / \n",
    "        #\n",
    "        \n",
    "        if result.success:\n",
    "            good +=1         \n",
    "        else:\n",
    "            bad += 1\n",
    "            \n",
    "        # 无论好坏，总是记录拟合结果, \n",
    "        with open(\"fitted_\" + random_name(), \"wb\") as FILE:\n",
    "            pickle.dump({\"result\":result, \n",
    "                         \"fitted_model\":fitted_model, \n",
    "                         \"lcs\":lcs, \n",
    "                         \"true_values\":[sn_z, sn_t0, sn_x0, sn_x1, sn_c]}, file=FILE)\n",
    "        #\n",
    "        # ---------------------------------------------\n",
    "\n",
    "#         print(f\"{sn_z}, {sn_t0}, {sn_x0}, {sn_x1}, {sn_c}\")\n",
    "#         print(result)\n",
    "    except (DataQualityError, RuntimeError):\n",
    "#         print(\"未观测到\")\n",
    "        bad +=1\n",
    "    \n",
    "print(f\"模拟的sn总数={good+bad}, 其中观测到{good}, 为观测到{bad}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "db2511d8-30b4-42cf-94fa-3bfedc45c05c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "模拟的sn总数=4675, 其中观测到2587, 未观测到309, 不满足信噪比条件的1779\n"
     ]
    }
   ],
   "source": [
    "print(f\"模拟的sn总数={good+bad+bad_input}, 其中观测到{good}, 未观测到{bad}, 不满足信噪比条件的{bad_input}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02691e24-0a9b-486e-a58f-e9d7d209aea7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  },
  "toc-showcode": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
